home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / fortran / linpkdrv.zip / SG.FOR < prev    next >
Text File  |  1984-01-05  |  18KB  |  616 lines
  1. C     MAIN PROGRAM
  2.       INTEGER LUNIT
  3. C     ALLOW 5000 UNDERFLOWS.
  4. C     CALL TRAPS(0,0,5001,0,0)
  5. C
  6. C     OUTPUT UNIT NUMBER
  7. C
  8.       LUNIT = 6
  9. C
  10.       CALL SGETS(LUNIT)
  11. C
  12.       STOP
  13.       END
  14.       SUBROUTINE SGETS(LUNIT)
  15. C     LUNIT IS THE OUTPUT UNIT NUMBER
  16. C
  17. C     TESTS
  18. C        SGECO,SGEFA,SGESL,SGEDI,SGBCO,SGBFA,SGBSL,SGBDI
  19. C
  20. C     LINPACK. THIS VERSION DATED 08/14/78 .
  21. C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
  22. C
  23. C     SUBROUTINES AND FUNCTIONS
  24. C
  25. C     LINPACK SGECO,SGESL,SGEDI,SGBCO,SGBSL,SGBDI
  26. C     EXTERNAL SGEXX,SMACH
  27. C     BLAS SAXPY,SDOT,SSCAL,SASUM
  28. C     FORTRAN ABS,AMAX1,FLOAT,MAX0,MIN0
  29. C
  30. C     INTERNAL VARIABLES
  31. C
  32.       REAL A(15,15),AB(43,15),AINV(15,15),ASAVE(15,15)
  33.       REAL B(15),BT(15),SDOT,DET(2),DETB(2)
  34.       REAL X(15),XB(15),XEXACT(15),XT(15),XTB(15),T,Z(15)
  35.       REAL AINORM,ANORM,SMACH,COND,COND1,EN,ENORM,EPS
  36.       REAL ETNORM,FNI,FNORM,ONEPX,RCOND,RCONDB,RNORM
  37.       REAL RTNORM,Q(8),QS(8),SASUM,XNORM,XTNORM
  38.       INTEGER I,IPVT(15),IPVTB(15),IQ(8),I1,I2,J
  39.       INTEGER K,KASE,KB,KBFAIL,KOUNT,KP1,KSING,KSUSP(8)
  40.       INTEGER L,LDA,LDAB,LUNIT,M,ML,MU,N,NM1,NPRINT
  41.       LOGICAL KBF
  42. C
  43.       LDA = 15
  44.       LDAB = 43
  45. C
  46. C     WRITE MATRIX AND SOLUTIONS IF  N .LE. NPRINT
  47. C
  48.       NPRINT = 3
  49. C
  50.       WRITE (LUNIT,460)
  51.       WRITE (LUNIT,880)
  52. C
  53.       DO 10 I = 1, 8
  54.          KSUSP(I) = 0
  55.    10 CONTINUE
  56.       KSING = 0
  57.       KBFAIL = 0
  58. C
  59. C     SET EPS TO ROUNDING UNIT
  60. C
  61.       EPS = SMACH(1)
  62.       WRITE (LUNIT,470) EPS
  63.       WRITE (LUNIT,450)
  64. C
  65. C     START MAIN LOOP
  66. C
  67.       KASE = 1
  68.    20 CONTINUE
  69. C
  70. C        GENERATE TEST MATRIX
  71. C
  72.          CALL SGEXX(A,LDA,N,KASE,LUNIT)
  73. C
  74. C        N = 0 SIGNALS NO MORE TEST MATRICES
  75. C
  76. C     ...EXIT
  77.          IF (N .LE. 0) GO TO 440
  78.          ANORM = 0.0E0
  79.          DO 30 J = 1, N
  80.             ANORM = AMAX1(ANORM,SASUM(N,A(1,J),1))
  81.    30    CONTINUE
  82.          WRITE (LUNIT,650) ANORM
  83. C
  84.          IF (N .GT. NPRINT) GO TO 50
  85.             WRITE (LUNIT,450)
  86.             DO 40 I = 1, N
  87.                WRITE (LUNIT,700) (A(I,J), J = 1, N)
  88.    40       CONTINUE
  89.             WRITE (LUNIT,450)
  90.    50    CONTINUE
  91. C
  92. C        GENERATE EXACT SOLUTION
  93. C
  94.          XEXACT(1) = 1.0E0
  95.          IF (N .GE. 2) XEXACT(2) = 0.0E0
  96.          IF (N .LE. 2) GO TO 70
  97.             DO 60 I = 3, N
  98.                XEXACT(I) = -XEXACT(I-2)
  99.    60       CONTINUE
  100.    70    CONTINUE
  101. C
  102. C        SAVE MATRIX AND GENERATE R.H.S.
  103. C
  104.          DO 90 I = 1, N
  105.             B(I) = 0.0E0
  106.             BT(I) = 0.0E0
  107.             DO 80 J = 1, N
  108.                ASAVE(I,J) = A(I,J)
  109.                B(I) = B(I) + A(I,J)*XEXACT(J)
  110.                BT(I) = BT(I) + A(J,I)*XEXACT(J)
  111.    80       CONTINUE
  112.             X(I) = B(I)
  113.             XT(I) = BT(I)
  114.             XB(I) = X(I)
  115.             XTB(I) = XT(I)
  116.    90    CONTINUE
  117. C
  118. C        FACTOR AND ESTIMATE CONDITION
  119. C
  120.          CALL SGECO(A,LDA,N,IPVT,RCOND,Z)
  121. C
  122. C        OUTPUT NULL VECTOR IF N .LE. NPRINT
  123. C
  124.          IF (N .GT. NPRINT) GO TO 110
  125.             WRITE (LUNIT,720)
  126.             DO 100 I = 1, N
  127.                WRITE (LUNIT,730) Z(I)
  128.   100       CONTINUE
  129.             WRITE (LUNIT,450)
  130.   110    CONTINUE
  131. C
  132. C        FACTOR BAND FORM AND COMPARE
  133. C
  134.          KBF = .FALSE.
  135.          ML = 0
  136.          MU = 0
  137.          DO 140 J = 1, N
  138.             DO 130 I = 1, N
  139.                IF (ASAVE(I,J) .EQ. 0.0E0) GO TO 120
  140.                   IF (I .LT. J) MU = MAX0(MU,J-I)
  141.                   IF (I .GT. J) ML = MAX0(ML,I-J)
  142.   120          CONTINUE
  143.   130       CONTINUE
  144.   140    CONTINUE
  145.          WRITE (LUNIT,790) ML,MU
  146.          IF (2*ML + MU + 1 .LE. LDAB) GO TO 150
  147.             WRITE (LUNIT,680)
  148.          GO TO 430
  149.   150    CONTINUE
  150.             M = ML + MU + 1
  151.             DO 170 J = 1, N
  152.                I1 = MAX0(1,J-MU)
  153.                I2 = MIN0(N,J+ML)
  154.                DO 160 I = I1, I2
  155.                   K = I - J + M
  156.                   AB(K,J) = ASAVE(I,J)
  157.   160          CONTINUE
  158.   170       CONTINUE
  159. C
  160.             CALL SGBCO(AB,LDAB,N,ML,MU,IPVTB,RCONDB,Z)
  161. C
  162.             IF (RCONDB .EQ. RCOND) GO TO 180
  163.                WRITE (LUNIT,780)
  164.                WRITE (LUNIT,820) RCOND,RCONDB
  165.                KBF = .TRUE.
  166.   180       CONTINUE
  167.             KOUNT = 0
  168.             DO 190 J = 1, N
  169.                IF (AB(M,J) .NE. A(J,J)) KOUNT = KOUNT + 1
  170.                IF (IPVTB(J) .NE. IPVT(J)) KOUNT = KOUNT + 1
  171.   190       CONTINUE
  172.             IF (KOUNT .EQ. 0) GO TO 200
  173.                WRITE (LUNIT,780)
  174.                WRITE (LUNIT,830) KOUNT
  175.                KBF = .TRUE.
  176.   200       CONTINUE
  177. C
  178. C           TEST FOR SINGULARITY
  179. C
  180.             IF (RCOND .GT. 0.0E0) GO TO 210
  181.                WRITE (LUNIT,710) RCOND
  182.                WRITE (LUNIT,480)
  183.                KSING = KSING + 1
  184.             GO TO 420
  185.   210       CONTINUE
  186.                COND = 1.0E0/RCOND
  187.                WRITE (LUNIT,500) COND
  188.                ONEPX = 1.0E0 + RCOND
  189.                IF (ONEPX .EQ. 1.0E0) WRITE (LUNIT,490)
  190. C
  191. C              COMPUTE INVERSE, DETERMINANT AND COND1 = TRUE CONDITION
  192. C
  193.                DO 230 J = 1, N
  194.                   DO 220 I = 1, N
  195.                      AINV(I,J) = A(I,J)
  196.   220             CONTINUE
  197.   230          CONTINUE
  198.                CALL SGEDI(AINV,LDA,N,IPVT,DET,Z,11)
  199.                AINORM = 0.0E0
  200.                DO 240 J = 1, N
  201.                   AINORM = AMAX1(AINORM,SASUM(N,AINV(1,J),1))
  202.   240          CONTINUE
  203.                COND1 = ANORM*AINORM
  204.                WRITE (LUNIT,510) COND1
  205.                WRITE (LUNIT,750) DET(1)
  206.                WRITE (LUNIT,760) DET(2)
  207. C
  208. C              SOLVE  A*X = B  AND  TRANS(A)*XT = BT
  209. C
  210.                CALL SGESL(A,LDA,N,IPVT,X,0)
  211.                CALL SGESL(A,LDA,N,IPVT,XT,1)
  212. C
  213.                IF (N .GT. NPRINT) GO TO 270
  214.                   WRITE (LUNIT,520)
  215.                   DO 250 I = 1, N
  216.                      WRITE (LUNIT,740) X(I)
  217.   250             CONTINUE
  218.                   WRITE (LUNIT,530)
  219.                   DO 260 I = 1, N
  220.                      WRITE (LUNIT,740) XT(I)
  221.   260             CONTINUE
  222.                   WRITE (LUNIT,450)
  223.   270          CONTINUE
  224. C
  225. C              MORE BAND COMPARE
  226. C
  227.                CALL SGBSL(AB,LDAB,N,ML,MU,IPVTB,XB,0)
  228.                CALL SGBSL(AB,LDAB,N,ML,MU,IPVTB,XTB,1)
  229.                KOUNT = 0
  230.                DO 280 I = 1, N
  231.                   IF (XB(I) .NE. X(I)) KOUNT = KOUNT + 1
  232.                   IF (XTB(I) .NE. XT(I)) KOUNT = KOUNT + 1
  233.   280          CONTINUE
  234.                IF (KOUNT .EQ. 0) GO TO 290
  235.                   WRITE (LUNIT,780)
  236.                   WRITE (LUNIT,840) KOUNT
  237.                   KBF = .TRUE.
  238.   290          CONTINUE
  239.                CALL SGBDI(AB,LDAB,N,ML,MU,IPVTB,DETB)
  240.                IF (DETB(1) .EQ. DET(1) .AND. DETB(2) .EQ. DET(2))
  241.      *            GO TO 300
  242.                   WRITE (LUNIT,780)
  243.                   WRITE (LUNIT,850) DETB
  244.                   KBF = .TRUE.
  245.   300          CONTINUE
  246. C
  247. C              RECONSTRUCT  A  FROM TRIANGULAR FACTORS , L AND U
  248. C
  249.                NM1 = N - 1
  250.                IF (NM1 .LT. 1) GO TO 330
  251.                DO 320 KB = 1, NM1
  252.                   K = N - KB
  253.                   KP1 = K + 1
  254.                   L = IPVT(K)
  255.                   DO 310 J = KP1, N
  256.                      T = -A(K,J)
  257.                      CALL SAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1)
  258.                      T = A(L,J)
  259.                      A(L,J) = A(K,J)
  260.                      A(K,J) = T
  261.   310             CONTINUE
  262.                   T = -A(K,K)
  263.                   CALL SSCAL(N-K,T,A(K+1,K),1)
  264.                   T = A(L,K)
  265.                   A(L,K) = A(K,K)
  266.                   A(K,K) = T
  267.   320          CONTINUE
  268.   330          CONTINUE
  269. C
  270. C              COMPUTE ERRORS AND RESIDUALS
  271. C                 E  =  X - XEXACT
  272. C                 ET =  XT - XEXACT
  273. C                 R  =  B - A*X
  274. C                 RT =  BT - A*XT
  275. C                 F  =  A - L*U
  276. C                 AI =  A*INV(A) - I
  277. C
  278.                XNORM = SASUM(N,X,1)
  279.                XTNORM = SASUM(N,XT,1)
  280.                ENORM = 0.0E0
  281.                ETNORM = 0.0E0
  282.                FNO